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1 . SUMMARY 

A method has been developed which calculates two- 

dimensional, transonic, viscous flow in ducts. The finite- 
volume, time-marching formulation is used to obtain steady 
flow solutions of the Reynolds-averaged form of the Navier 
Stokes equations. The entire calculation is performed in the 
physical domain. 

The features of the current method can be summarized as 
follows. Control volumes are chosen so that smoothing of 
flow properties, typically required for stability, is not 
needed. Different time steps are used in the different 

governing equations. A new pressure interpolation scheme is 
introduced which improves the shock capturing ability of the 

method. A multi-volume method for pressure changes in the 

boundary layer allows calculations which use very long’ and 
thin control volumes (length/height - 1000). The method is 
then compared here with two test cases. Essentially incom- 
pressible turbulent boundary layer flow in an adverse pres- 
sure gradient is calculated and the computed distributions of 
mean velocity and shear stress are in good agreement with the 
measurements. Transonic viscous flow in a converging diver- 
ging nozzle is calculated; the Mach number upstream of the 
shock is approximately 1.25. The agreement between the 
calculated and measured shock strength and total pressure 
losses is good. 

2. INTRODUCTION 

The finite volume method has been used extensively to 
solve the Euler equations for transonic flow including flow 
at high Mach numbers. In internal aerodynamics, McDonald [1] 
was the first investigator to use the time marching finite 
volume method. Denton [2] extended McDonald's finite-volume 
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method to three dimensions. Versions of Denton's method have 
been used in Inviscid-viscous interaction programs for 
turbomachinery calculations [3-5], 

The scope of the present work was to extend a finite 
volume method like that of Denton's to be able to calculate 
laminar or turbulent flow in ducts. The new method has the 

capability to calculate subsonic as well as transonic flow. 

3. GOVERNING EQUATIONS 

The unsteady form of the continuity equation, the x- 
momentum equation, and the y-momentum equation, in integral 
form, are used to obtain a steady-state solution for flow 
through 2-dimensional ducts. The ideal gas equation of 
state, the assumption of constant total temperature, and a 
Prandtl mixing length turbulence model complete the governing 
equations needed to solve for the unknown variables p, u, v, 
P, u , and T. 

For a finite control volume where we can assign one 
value of density to the control volume, and for a finite time 
step, St, continuity states that, 

P n+1 - p" - Sp » -[// pu • dA] (1) 


where the integral is evaluated explicitly at the current 
time step, n. In arriving at an expression which relates the 
pressure change directly to the continuity error, we will 
assume that changes in temperature are small in comparison to 
other changes for one time step. Thus, we can relate changes 
In pressure to changes in density through the ideal gas 
equation of state, 

p n+l - p n , 5P » -r T [JJ P u . <JA] (2) 

For the method introduced in the current work, a non-conser- 
vative form of the unsteady momentum equation is used. The 
non-conservative form is used because it allows the use of 
different time steps for the continuity and momentum equa- 
tions. The differences between the non-conservative and 
conservative forms of the unsteady momentum equations are 
associated with the unsteady and convective terms. Speci- 
fically, we note that 

3(pu) 3u 

" + V • pu u » p .jZ + pu • Vu (3) 
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and the right hand side of Eq. (3) can be rewritten as 


3u 

p + pu • Vu 


3u 

p 7t + V 


pu u - u(V 


pu) 


( 4 ) 


When the right hand side of Eq. (4) is combined with the 
pressure and viscous terms, the momentum equation in integral 
form becomes 


(u) n+1 - (u) n ■ 6(u) ■ [-// pu u • dA + u // pu • dA 

“// P * dA + // (uVu + yVu T • dA)] (5) 

To maintain stability, the properties must be updated in the 
proper sequence. In the current method, the sequence is 

1. update the pressure from the continuity equation; 

2. update the velocities from the momentum equation using 
the new pressure and old velocities and old density; 

3. update the density from the ideal gas equation of state; 

4. update the temperature from constant total temperature. 

4. CONTROL VOLUMES 

A new control volume has been introduced for this 
method. To eliminate the need for smoothing of flow proper- 
ties, there must be as many control volumes across the duct 
as there are nodes where these variables are calculated. We 
need as many equations as unknowns. The control volumes also 
need to be located so that errors in continuity and momentum 
can correctly influence the changes in pressure or density 
and velocity without smoothing. The current control volume 
accomplishes this and is shown in Fig. 1. When calculating 
the flux through a streamwise face of an element, the value 
of the flow properties at the node on that face are used. 
Linear interpolation is used to obtain the flux on the cross- 
stream face. 



Fig. 1 New Control Volumes 
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5. DISTRIBUTION OF PROPERTIES 


The properties at node points are changed in the flow 
field after each time step because the continuity and 
momentum equations are not satisfied for a given control 
volume. A decision must be made about which node, either 
upstream or downstream, these changes should be allocated 
to. The criterion used in determining where changes in 
properties should be sent is that these distributions result 
In reduced errors in continuity and momentum. The current 
method uses the following allocation procedure: 

1. The pressure is updated through the continuity equation 
and the pressure change is sent to the upstream node; 

2. The u and v velocities are updated through the momentum 
equations and the changes are sent to the downstream 
node; 

3. The density is updated through the ideal gas equation of 
state using an interpolated pressure. 

6. PRESSURE INTERPOLATION PROCEDURE 

As part of the updating procedure used by Denton [5], an 
effective pressure is used in the momentum equations rather 
than the true thermodynamic pressure determined from the 
Ideal gas equation of state. This effective pressure is 
needed because If the true pressure is used in the momentum 
equations the solution may not converge. In the current 
method, the density used In the continuity and momentum 
equations is an effective density which may be different than 
the density obtained using the ideal gas equation of state. 
This effective density is used to satisfy stability 
requirements. 

Starting with a generalized pressure interpolation 
equation for the effective density 


(P - P ) 

P I+1 * t P I + a 0 (P I+l " P I ) + a l 2 + 

( P I+1 “ P I— 2^ i 1 

a 2 3 J ’ 

Mach number limitations were sought for a Q , aj, and a 2 
that 


( 6 ) 

such 


a l + a 2 + a 3 = 1 (?) 

which assures second order accurate solutions. A set of 
equations for ag, aj, and a 2 was chosen which satisfies two 
stability criteria [6], The equations are 
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for M < 2 


; a x - 


a 0 " 



(\" 1 ) 

vr 


1 - aQ ; &2 * 0 : 


for M > 2 


0 ; a. 


4/M^ 


1 - a, 


( 8 ) 


These Mach number dependent formulations for ag, a^, and &2 
are shown in Fig. 2. 



HACK HUMBER 

Fig. 2 Mach Number Dependent Values for 
Coefficients ag, a^, and &2 


7. TIME STEPS 

A unique feature of this method is the use of different 
time steps for the continuity and momentum equations. 
Previous workers who have used explicit time marching methods 
have used the CFL condition as a basis for determining 
allowable time steps which maintain stability. The same time 
step is used for both the continuity and momentum 

equations. In the current method, the expressions that are 
used to determine the allowable time steps are; for the 
momentum equations 


u 

4. 

v eff 

+ 1 2u 1 

Tx 


«5y 

'p(Sy) 2 ' 


(9) 


and for continuity. 


St < 
c 



1 



U J 

V 

mx| + 



( 10 ) 


where St is the momentum time step, 5 1 is the continuity 
time step and v e £f is an effective y-component of velocity. 
The advantage of using different time steps is that, for low 
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velocity regions of the flow, the allowable momentum time 
step can be significantly larger than that allowed by the CFL 
condition. These larger time steps allow the boundary layer 
profiles to change more rapidly and enhance the convergence 
rate significantly compared with a method which uses the CFL 
condition. 

8. BOUNDARY CONDITIONS 


For viscous flow, at the upstream boundary, the total 
temperature, freestream total pressure, inlet boundary layer 
velocity profile, and flow angle are specified. Along the 
downstream. boundary the static pressure is specified. Pres- 
sures along the solid boundaries are determined from linear 
extrapolation. For viscous flow, the values of the x- 
component and y-coraponent of velocity are set equal to zero 
at solid walls. 


9. TURBULENCE MODEL 


A Prandtl mixing length model is used to model the 
turbulent stresses. The model is 


U eff 


+ M , 



"du" 

‘3y 


L is smaller of 0.08 times the width of boundary layer 
or 0.41 times the distance to the wall 


Van Driest Correction 

L - 0.41 "y"(l - exp[- "y" /pT/26 u £ ]) 
Near Wall Correction ^eff * + h t ) 


10. MULTI-VOLUME METHOD FOR PRESSURE CHANGES 

Control volumes are grouped in the boundary layer to 

form a larger global control volume. The continuity error is 

calculated for this global control volume and changes in 

pressure are assigned equally to each of the upstream nodes 
for each control volume making up the global control 
volume. Then the global control volume is made successively 
smaller towards the wall. This is shown schematically in 

Fig. 3. The entire pressure change for one Iteration at each 
node within the multi-volume region is determined by adding 
together all the pressure changes assigned to that node. 

The multi-volume method propagates pressure changes 
rapidly through the boundary layer and minimizes transverse 
pressure gradients in the intermediate solution. The above 
changes allow the calculation of boundary layer flows where 
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the control volumes near the wall can have aspect ratios 
(length/height) over 1000. 



VOL. 1 VOL. 2 V0L - 3 VOL. 4 

Fig. 3 Multi-Volume Method for Pressure Changes 
in the Boundary Layer 

11. TRANSVERSE UPWIND DIFFERENCING 

When the control volumes become long and thin near the 
wall of the duct, the fluxes through the top and bottom faces 
of the control volume become more significant in comparison 
to the fluxes through the streamwlse faces. To strengthen 
the diagonal dominance of the coefficient matrix, the 
momentum fluxes through the transverse faces may be calcu- 
lated using interpolated velocities upstream in the 
transverse direction rather than the actual interpolated 
values. The Interpolation functions and the derivation of 
the functions is discussed in more detail in Ref. 6. 

12. SAMUEL AND JOUBERT INCOMPRESSIBLE TURBULENT BOUNDARY 
LAYER 

Incompressible turbulent boundary layer flow in a 
diverging duct was calculated for test case 0141 of the 

Stanford Conference [7]. The grid used in the present 
calculations is shown in Fig. 4. The inlet velocity is 26 
m/s. 



Fig. 4 Geometry and Grid for Samuel and Joubert 

Figure 5a shows a comparison of the calculated skin 
friction coefficient with the experimental results and with 
the results from the Moore cascade flow program. The 
agreement is excellent. A comparison of the calculated 
turbulent shear stress distribution, uv, with the experi- 
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mental results is shown in Fig. 5b. The agreement is good. 
Figure 5c shows good agreement also between the calculated 
and measured velocity profiles at two locations in the duct. 
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Fig, 5 Results for Samuel and Joubert 
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13. MDRL DIFFUSER CALCULATIONS 


The diffuser geometry (Model G) Is shown In Figure 6a 
[8,9]. Figure 6a also shows the computational grid used 
which has 87 grid points in the axial direction and 20 points 
across the flow. The Inlet boundary layer thicknesses were 
specified as 9% and 4.5% of the inlet diffuser height for the 
curved and flat wall boundary layers, respectively. For this 
calculation, the ratio of exit static pressure to the inlet 
total pressure was 0.826. In the experiment, this test point 
results in transonic flow in the diverging portion of the 
duct with a Mach number of approximately 1.235 upstream of a 
nearly normal shock, and the flow remained fully-attached 
throughout the diffuser at this test condition. 

A contour plot of static pressure is shown In Fig. 6b. 
The shock can be seen in the diverging portion of the duct. 
The shock Is well defined as illustrated by the high 
clustering of contours at the shock. Figure 6c shows a Mach 
number contour plot for the calculations. The calculated and 
measured curved wall static pressures are compared in Fig. 
7. The shock is well defined and no overshoot occurs in the 
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b) static pressure contours 
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c) Mach number contours 


Fig. 6 Geometry and Contours for MDRL Diffuser 

Measured shock locations on the curved wall and in the 
middle of the duct are plotted in Fig. 8 as a function of 
shock Mach number, M au , determined from the minimum wall 
static pressures on the curved wall. The minimum wall static 
pressure in the calculation is located at x/h “ 1.5; this is 
taken to be the location of the shock. The Mach number 
upstream of the shock was determined to be 1.256 from the 
calculated total pressure ratio across the shock in the 
frees tream. This result is plotted in Fig. 8 and it agrees 
well with the measured shock location. Comparisons of 
calculated and measured velocity profiles (see Ref. 9) at two 
axial locations along the duct are shown in Fig. 9. The 
agreement is good. The mass averaged total pressure at the 
diffuser exit divided by the inlet freestream total pressure 
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Fig. 7 Curved Wall Static Pressures for MDRL Diffuser 



CALCULATED 


Fig. 8 Comparison of Computed and Measured 
Shock Position in MDRL Diffuser 



+ CALCULATION 
A EXPERIMENT 


«.j i.i e.« m 



+ CALCULATION 
A EXPERIMENT 


Fig. 9 Velocity Profiles at x/h= 4.03 and 8.2 
in MDRL Diffuser 



is calculated from the numerical results to be 0.9615. This 
compares well with the measured value of 0.965, obtained from 
the data of M. Sajben and T. J. Bogar, midway between the 
diffuser side walls. 

The total CPU time for the MDRL diffuser calculations 
was approximately 35 minutes on an IBM 3031. 

14 . CONCLUSIONS 

An explicit finite volume time marching method has been 
extended to allow the calculation of laminar and turbulent 
flow in ducts. Both subsonic and supersonic flow can be 
calculated with the method. Incompressible turbulent 
boundary layer flow in an adverse pressure gradient was 
calculated. The agreement between the calculated and 
measured skin friction coefficient, turbulent shear stress 
distribution, and mean velocity profiles was good. Transonic 
viscous flow through a converging diverging nozzle was 
calculated. The computed and measured velocity profiles were 
in good agreement especially near the exit of the nozzle. 
The computed and measured shock locations were compared and 
were found to be in good agreement. Viscous and shock losses 
in the diffuser were well modelled. 
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